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TECHNICAL PAPER 


A STOCHASTIC MODEL FOR PHOTON NOISE INDUCED BY CHARGED PARTICLES 
IN MULTIPLIER PHOTOTUBES OF THE SPACE TELESCOPE 
FINE GUIDANCE SENSORS 


I. INTRODUCTION 


A. General 

The pioneering observatory Space Telescope (ST) is scheduled for operation in Earth-orbit 
beginning in mid-1986. 

For the ST to perform its mission, very accurate line-of-sight stability must be maintained for 
periods up to 10 hr. To maintain this stability, the attitude information is provided by two on-board 
Fine Guidance Sensors (FGS) (a third FGSS is either on stand-by or it is performing astrometry). The 
FGS track guide-stars through the main telescope. Each FGS has two multiplier phototube (PMT) pairs 
(Fig. 1), one pair for each axis. Pointing stability is maintained by detecting deviations from the viewing 
object using the difference in guide star photon counts in one or more of the PMT pairs. These photons 
are detected by a photocathode on the inner side of the PMT glass window. Because of ST exposure to 
charged particles in its space environment, concern arose over the possibility that incoming particles 
would strike the glass window, creating a flash of Cerenkov radiation 1 and secondary fluorescence as the 
particle bores its way through the glass. 
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Figure 1. Multiplier phototube cross section. 

1. Cerenkov radiation. The radiation from a charged particle whose velocity is greater than the phase 
velocity that an electromagnetic wave would have if it were propagating in the medium. The particle 
will continue to lose energy by radiation until its velocity is less than this phase velocity. This phenom- 
enon is analogous to the generation of a shock wave when an object is traveling faster than the velocity 
of sound in the medium. A bow wave is set up that radiates energy into the medium and slows down 
the object. The angle that the cone of radiation makes with the direction of motion of the particle can 
be used to measure the velocity of the particle. Named after the Russian physicist Pavel A. Cerenkov, 
born 1904; Nobel prize 1958 [1], 




Because the photocathode cannot distinguish the source of the photons, those photons generated 
by the charged particle (which then strike the photocathode) will cause an increased photon count in the 
struck PMT. This would appear as an increase in apparent star brightness and would produce a large 
differential count between pairs of PMTs resulting in a ST attitude disturbance, possibly large enough 
to lose “lock” on the desired guide star. The Cerenkov flash is very short lived (less than 70 ps) and 
causes only one photon count with an associated pulse amplitude proportional to the number of 
Cerenkov photons striking the photocathode because the time resolution of the PMT pulse amplitude 
discriminator (PAD) is approximately 75 ns. However, the duration of the fluorescence follows an 
exponential distribution, with a time constant of several 100 jus and every fluorescence photon is 
counted. This gave rise to the idea to detect the Cerenkov flash by an additional much higher PAD 
threshold and then blank the photon count for the next millisecond, eliminating the effect of the 
fluorescence photons. While others began to study the problem from an experimental point of view, we 
decided to use a simulation to derive the bivariate probability distribution for the number of Cerenkov 
counts and for the number of fluorescence counts registered by the photocathode for a single charged 
particle strike. 

Because of the complexity of the physical setting and apparent random behavior of particle 
strikes, it was decided that a stochastic model for photon noise induced by charged particles be 
developed based on a computer simulation of the physical scenario. The resultant bivariate distribution 
of Cerenkov radiation and fluorescence would show the relationship between these two and determine 
whether the Cerenkov flash can be used to give reliable warning of the subsequent fluorescence. Secondly, 
these results may then be introduced into the ST attitude-control model to determine the effects of 
charged particles on ST pointing accuracy and to suggest possible remedies if required. 

A comparison of the model’s predictions with experimental results is also discussed, as well as a 
sensitivity study of some of the model’s parameters. 


B. Pointing Capability 

The ST has three sources of information for its orientation: the Rate Gyroscopes, the Star 

Trackers, and the Fine Guidance Sensors. Since the primary goal is the prompt and accurate acquisition 
of preselected guide stars, all three sources of information are used when ST is pointed to a new object. 
The slewing rate is determined by reaction wheels and is specified to allow movement to objects 90 deg 
away in no more than 20 min, including angular acceleration and deceleration, with the reference frame 
provided by the Rate Gyroscopes. The Star trackers then use bright stars to determine the pointing to 
an accuracy of a few arc minutes, which is sufficient to place the guide stars in the field of the Fine 
Guidance System. Following an acquisition procedure to select the correct guide star, the guide star 
images are put through a 3 arc sec x 3 arc sec aperture and fed to an interferometric device which gives 
(through the use of a pair of PMTs per axis), the fine error signal that results in the overall guiding 
stability of 0.007 arc sec rms. The positioning of the Fine Guidance Sensors can be set to an accuracy 
of 0.001 arc sec. 


C. Charged Particle 2 Environment 

The nominal orbit for the ST is at an inclination of 28.8 deg and an altitude of 590 km (320 
n.mi.). For this orbit, radiation-sensitive equipment will be subjected to a substantial flux of high energy 
charged particles and their induced secondary particle showers. 

2. In this report, the term charged particle will apply to both cosmic ray nuclei and those particles 
trapped in the Van Allen Belt by the Earth’s magnetic field. 


2 



There are three possible sources of radiation that can degrade the ST pointing performance: geo- 
magnetically trapped electrons, protons, alpha particles (helium nuclei), and heavier nuclei, solar event 
protons, and galactic cosmic ray particles. Because the Earth’s magnetic field at 30 deg latitude will 
eliminate protons below 4 BeV, solar flare protons should not be important. The dose rate of charged 
particles will be reduced by the geomagnetic field and the vehicle structure will shield against electrons. 
According to Stassinopoulos [2], the dose rate caused by particles for the nominal orbit should be about 
0.7 rads/day (a rad is the deposition of 100 ergs of energy per gram of material receiving the dose). 

When a relativistic charged particle passes through the multiplier phototube window, its speed is 
greater than that of light in the glass. This produces a cascade of conically directed photons, the resultant 
flash of Cerenkov radiation lasting less than 70 ps. Because the time resolution of the PMT electronics is 
75 ns, the Cerenkov radiation is registered only as a single pulse of amplitude proportional to the number 
of Cerenkov photons striking the photocathode. During the next several hundred jus, fluorescence photons 
are emitted and those which strike the photocathode are individually registered as photon noise. 

Since the photocathode is incapable of distinguishing between these photons and the incoming 
guide star photons, the biased 25 ms (millisecond) interval (cycle length of ST computer) count will be 
interpreted by the ST guidance system as an attitude offset from the guide star which will result in 
pointing disturbances. 


D. Modeling Approach 

The rate at which charged particles would strike the PMT window had been estimated from 
satellite data by others. Hence, it is sufficient to determine the resultant radiation effects of each strike. 
Consequently, the distribution of Cerenkov radiation and fluorescence is to be simulated from the physical 
scenario, conditional on having a charged particle strike the PMT window. Given that such a charged 
particle had struck the PMT window, its entry point was simulated by assuming that the likelihood of 
it striking the top, side wall, or bottom, would be proportional to that region’s surface area, where, for 
geometrical convenience, the PMT window was idealized as a right circular cylinder. The outer window 
surface was considered the “top” and the inner surface the “bottom”. The surface of entry was chosen 
by Monte Carlo methods and then a uniformly randomly distributed point on that surface was designated 
as the entry point. 

Because direction of entry was assumed to have no preference, a uniformly distributed random 
point was selected on the surface of a hemisphere imagined to be placed over the entry point. The two 
points were then connected, giving the direction cosines of the particle’s path as it enters the PMT 
window. The particle was followed to its exit point and the path length determined. Instantaneously 
Cerenkov photons are emitted, their number being a random variable assumed to follow a Poisson dis- 
tribution with a mean proportional to the particle’s path length through the glass and origin of departure 
of the individual photons being uniformly randomly distributed along the path. The direction of travel 
was randomly selected down the surface of a cone whose half-angle is determined by the particle’s kinetic 
energy, with apex at the photon departure point and base towards the charged particle’s exit. The distri- 
bution of the kinetic energy of charged particles trapped in the Van Allen Belt would also be estimated 
from available satellite data for the ST orbit, while a literature search revealed that the cosmic ray energy 
distribution followed a power law [3,4]. 

The path of each of these photons was then followed to determine its point of contact with the 
cylinder’s interior surface. Those which were intercepted by the photocathode on the bottom left the 
system and contributed to the pulse amplitude of the Cerenkov radiation. Photons striking other regions 
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of the cylinder’s surface either passed through the cylinder’s surface and left the system or were reflected 
back into the cylinder depending on the angle of incidence. Reflected photons eventually would leave the 
system through either the photocathode or cylinder walls, or fall into trapped orbits within the cylinder 
and thus would not contribute to the FGS photon count. 


II. MATHEMATICAL MODEL 


NASA scientists had estimated that the probability of more than one cosmic ray strike in a 25 
ms interval (the basic computational cycle of the ST computer) would be negligible. 

Thus, this model is concerned with simulating the discrete bivariate distribution for the number 
of photons striking the photocathode from the Cerenkov radiation and the secondary fluorescence, 
respectively, given the event of a single charged particle having struck the PMT window. 

The PMT window was idealized as a right circular cylinder with radius R and height H (Fig. 2). 
The effective photocathode was represented as a disk of radius r on the inner (bottom) side of the PMT 
window. Since every surface element of the PMT glass window is exposed to 2n steradians of space, it 
was assumed that the probability of a hit on any one of the cylinder’s distinct surfaces (top, bottom, or 
side) was proportional to the surface’s area. These areas are 


tb = * r2 

top or bottom area 

‘ Sd = 2ffH R 

side area 

tot = 2 * R ( R + H) 

total area 


with R being the glass window radius, and H the glass window height. (See Appendix F for values of model 
parameters; Fig. 2 not drawn to scale.) 


z 



Figure 2. PMT window geometry. 
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The probabilities are then: 


P t b = \bl\ct - R /« + H) 


(l) 


p sd * A sd/ A tot = « /(R + H) . 


( 2 ) 


The direction cosines, selected such that the direction of the particle’s path is uniformly distri- 
buted over 4w steradian, are given by 3 


[dp d 2 , d^]^ = [cos <j> cos 0, cos sin 0, sin0]T 


(3) 


where 


6 = sin' 1 [2(U 2 - 0.5)] 


elevation angle 


0 = 2ttU 2 


azimuth angle, 


and the Us are different random numbers uniformly distributed between zero and one. 
A decision is now made, whether the top/bottom, or side was hit: 


U 


< P 


tb 


> P 


tb 


yields a top/bottom hit 


results in a side hit. 


A. Top/Bottom Hit Entry 

The coordinates of the top/bottom entry point are determined by 


x = R U , 
y = RU , 

(the Us are different uniform random numbers) and 

3. See Appendices A, B, and C for the proof. 


(4) 


5 



I H for d 3 < 0 

0 for d 3 > 0 . 


The requirement that points (x,y) be uniformly distributed over a disk of radius R is met by 
first generating points uniformly distributed over the first quadrant of edge R and then rejecting those 
which lie outside the circle of radius R. (Because of symmetry, only one quadrant is necessary. For 
graphical purposes, the quadrants were later randomly indexed.) 


B. Exit for Top/Bottom Entry 

The change in position from entry point to exit point is first calculated as if the top or bottom 
of the window were an infinite plane. Thus, 


"Ax' 


" d l 



Ay 

= 

d 2 

[H I(d 3 )-z]/d 3 

(5) 

Az 


_ d 3 . 




where I is the indicator function defined by 

I I for dg > 0 

0 for d 3 < 0 


For top or bottom entries, z is either H or 0, but the same equation can be used for starting points after 
reflections, where z can be between 0 and H. 

Notationally, any time a path or partial path is started at the window surface, the coordinates 
of its position are labeled as (x,y,z). Path endpoints are labeled (u,v,w). Ending can occur as an exit from 
the window, or with reflection, as the start of the next partial path. 

In order to determine the exit point, the change in position is added to the entry point 
coordinates: 


u 


X 


Ax" 

V 

= 

y 

+ sc 

Ay 

w 


z 

to to 


Az 
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A check is now made, whether this exit point (u,v,w) is within the window radius (scaling 
constant SC=1). If not, SC is reduced so that this condition is met by 


SC = (b + 



+ a(R^ - x^ - y^) )/a 


(7) 


where 


a = (Ax)^ + (Ay)^ 


and 


b = -(xAx + yAy) 


C. Side Hit Entry 

Since the window is symmetric, the only relevant parameter for a side hit is the height z above 
the bottom of the window (z=0) and it is determined to be uniformly distributed as 


z = H U 


( 8 ) 


The other coordinates are arbitrarily fixed at 


x = 0 and y = -R . 


However, for graphical purposes, we then rotate side entries through a uniformly distributed 
angle between 0 and 2n radians about the z-axis. 

As previously discussed, direction cosines of the charged particle’s path could be generated by 
connecting its entry point to a point uniformly distributed over the surface of a hemisphere centered 
over the entry point. However, mathematically, we found it more convenient to generate points distri- 
buted over the surface of a sphere and use the sign of the third direction cosine (d 2 ) to distinguish top 

and bottom particle strikes [equation (4)]. We set d 2 = |d 2 1 in order to direct side hits always into the 
window. 


D. Exit for Side Hit Entry 

The unsealed exit position is 


u = 2R dj d 2 /(l - dg^) 

v=(ud 2 /d 1 )-R (9) 

w = z + Az , 
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where 

Az = u d /d 

To make sure that the path is physically possible, the following scaling checks have to be made 

1) No top scaling is necessary if Az < (H - z) 

2) No bottom scaling is necessary if -z < Az. 

For top scaling we get 

SC = (H - z)/Az 
u = SC u 

( 10 ) 

v = SC(v + R) - R 
w = H . 

For bottom scaling we use 

SC = -z/Az 
u = SC u 

(ID 

v = SC (v + R) - R 
w = 0 

E. Length Determination 

Given the particle’s entry and exit points, the path length 

L = yj (x - u) 2 + (y - v) 2 + (z - w) 2 (12) 

is needed to determine the average number of Cerenkov and fluorescence photons. 
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F. Charged Particle Energy and Relative Speed 


We distinguish two types of charged particles according to their energy distribution: Cosmic ray 

nuclei and charged particles trapped in the Earth’s radiation belts. A power law (proportional to E"?) for 
the cosmic-ray energy spectrum is assumed, with the values of 7 appearing in Table 1, where Z is the 
charge number or atomic number [3,4]. 

TABLE 1. CHARGED PARTICLE SPECTRAL INDICES, REST MASSES 
AND RELATIVE ABUNDANCES 


z 

Particle 

Index 

Mass 

(kg) 

Relative Abundance 
Above 1 GeV/Nucleon 
(m^ . sr • s)*l 

1 

Proton 

2.70 

1.67E-27 

700 

2 

Alpha particle 

2.70 

6.6E-27 

140 

3 

Lithium 

3.35 

1 1.5E-27 

0.2 

4 

Beryllium 

3.35 

15.0E-27 

0.2 

5 

Boron 

3.35 

18.0E-27 

1.5 

6 

Carbon 

2.76 

20.0E-27 

4.9 

7 

Nitrogen 

3.07 

23.0E-27 

1.4 

8 

Oxygen 

2.67 

26.7E-27 

4.5 

9 

Fluorine 

3.01 

31.7E-27 

.18 

10 

Neon 

2.88 

33.1E-27 

.73 

11 

Sodium 

3.14 

38.4E-27 

.21 

12 

Magnesium 

2.60 

40.6E-27 

.70 

13 

Aluminum 

2.52 

45.1E-27 

.11 

14 

Silicon 

2.63 

46.9E-27 

.53 

26 

Iron 

2.45 

93.3E-27 

.08 


Thus, the following family of probability density functions for cosmic-ray energies per nucleon 
are obtained: 


f(E) = (7-I) E-r 


( 13 ) 


where 1 < E < 00 (BeV), which leads to 


e = u l i l ~y 
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as the desired random number generator. Because of the glass medium, nuclear interactions do not occur 
with significant frequency [5], Thus, the more complicated electromagnetic cascade process referred to 
as the “soft component” of cosmic rays by Janossy [6], in which photons give rise to energetic electrons 
which in turn give rise to photons, is not modeled. The electromagnetic cascade process eventually comes 
to an end because of the loss of energy by ionization. 

For protons trapped in the radiation belts, the flux differential for the ST orbit is shown in 
Figure 3 [7]. Because these protons are of much lower energy than the cosmic ray protons, the effect 
of shielding must be considered. The effective shielding for the ST was taken to be 7 mm of aluminum 
corresponding to the curve labeled by X=20 in Figure 3. A log-normal probability density function was 
found to be adequate, as suggested by References 8 and 9. It is of the form 


f(E) = 



“7 (In E - (x) 2 

e 2 * 2 


(14) 


where 


0 < E < oo 


The desired log-normal generator for Van Allen belt proton and alpha particle energy per nucleon is: 


E = e ix + ax 


(15) 


where x is a standard normal random number and n = 4.7 and a = 1. 



10 


Figure 3. Average trapped proton flux differential with respect to energy at the center 
of a spherical aluminum shell of the given shield thicknesses for the ST orbit. 



Thus, when modeling charged particle strikes in the Earth’s radiation belts, equation (15) is used 
to simulate the particle’s energy and equation (13) is used for cosmic ray nuclei strikes. Energies are then 
multiplied by Kj FBEV in order to convert to the SI system of units, where K JFBE y = 1.60219E-10 is 

the conversion factor from billion electronvolts to Joules. The charged particle’s speed relative to the 
speed of light c = 2.99792458E8 m/s can now be determined by 


P = 



- 1/(1 + E/m 0 c 2 ) 


(16) 


where the primary charged particle’s rest mass m Q is obtained from Table 1. 


G. Number of Induced Photons 

The Poisson mean for the number of induced Cerenkov photons is [10] 


N cr = (27rZ 2 /k) [1 - 1/((3N g ) 2 ] F(A) L 


(17) 


where Z is the charge number of the particle (atomic number); k = (47t//uq) e 2 /hc = 137 is the dimensionless 
fine structure constant, m 0 is the permeability coefficient, and h is Plank’s constant, N G =1.49 is 

the refractive index of 7056 Pyrex glass; and L is the path length (mm) computed by equation (12). 
The factor F is a function of the wavelength of light and is given by F(A) = 244 for 7056 Pyrex glass 
which was derived as follows [equation (2) of Reference 10] : 


F(\) = 1/Aj - 1/X 2 

for constant quantum efficiency. For varying efficiency e(A), F(A) is expressed as the integral 

^2 

F(A) = -0.5 j [e(X)/X 2 ] dX . (18) 


Figure 4 shows the efficiency of 7056 Pyrex glass as a function of the wavelength [11]. Equation (18) 
was numerically evaluated from Xj = 350 nm to X 2 = 800 nm, resulting in 


F(A) = 0.000 244 1/nm = 244 mm'* 
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Figure 4. Radiant energy spectral coverage 

The number of resultant photons comprising the fluorescence was assumed to follow a Poisson 
distribution with mean [12] 


n FL = K FL l ( dE / dx ) 


(19) 


where 


KpL = 1 6 for 4 tt steradian , [13] 


and dE/dx is determined by the restricted energy loss distribution (Appendix D). 

Hence, the simulated number of induced photons for both the Cerenkov pulse amplitude height 
and secondary fluorescence is determined by the inverse cumulative distribution method described in 
Appendix B using the Poisson cumulative distribution function with the above means of and Npp, 
respectively. 

The simulated number of Cerenkov photons n^p and fluorescence photons npp are then modeled 

as being uniformly randomly distributed along the charged particle’s path. Their coordinates are deter- 
mined by 
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(XpYpZj) = (x,y,z) + (Ax,Ay,Az) Uj , i = 1, n CR 

for Cerenkov photons and 

(Xj^zj) = (x,y,z) + (Ax,Ay,Az) Uj , i = 1, n FL 

for fluorescence, where the U’s are uniformly distributed random numbers and Ax, Ay, and Az are the 
components of the pathlength L. 

For Cerenkov radiation, the angle that the cone of luminescence makes with the direction of 
motion of the charged particle is determined by References 10 and 12 to be 

6 = cos' 1 (1/j 3N g ) (20) 

where |3 is determined from equation (16). The direction of departure of each of the n GR photons will 

make an angle 6 with the charged particle’s path. The three-dimensional effect of the luminescent cone 
is simulated by rotating each photon’s path through a random angle uniformly distributed over (0,27r) 
using the Euler transformation 


” d l~ 


cos 0:3 

-sin a 3 

o" 


cos «2 

0 

sin <*2 


cos 0:4 

sin 0:4 

o - 


ro - 

d 2 

= 

sin 0:3 

cos (X3 

0 


0 

1 

0 


-sin «4 

cos 0:4 

0 


sin oi] 

_ d 3 _ 


0 

0 

1 


-sin «2 

0 

cos <*2 


0 

0 

1 


cos cej 


( 21 ) 


with 

= cos' 1 (l/j3N G ) 

<*2 = cos' 1 (dj) 

«3 = tan' 1 (dj^) 

0:4 = 27rU (U is a uniform random number) . 


Fluorescence photons are emitted uniformly in all directions. Consequently, direction cosines 
are generated for each of the n FF photon paths by generating points uniformly distributed over the 

surface of a sphere and connecting each to the photon’s departure point. 
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H. Total Reflection 


The intensity of the reflected light R is [14] 

Nj cos - Nj cos 0 t 
Nj cos 0j + N t cos 0 t 

and 

Nj. cos 0j - Nj cos 0^. \ 2 
Nj cos 0 t + N t cos 0^ J 

where the intensity of the incident light is unity for perpendicularly polarized light (1) with respect to 
the plane of incidence as well as for the parallel polarized beam ( II) at the interface from a denser optical 
medium N| to a less dense and 0j is the angle of incidence, 0^ the angle of transmission, both with 

respect to the interface normal. Figure 5 shows and R|| as a function of 0 using Nj = Nq = 1.49 

and N t = 1 . From an angle of 42.4 deg to 90 deg total reflection is obtained. Also shown in Figure 5 is 
the average 





R = (Rj. + R||)/2 


(23) 


which is used in the model under the assumption that the light is not necessarily polarized. The general 
geometry of reflection is shown in Figure 6. Because individual photons in the model are dealt with, R 
is interpreted to be the probability of being reflected. If a standard uniform random number is greater 
than R, the photon is lost through transmission. Otherwise it is reflected and then traced to its next 

point of contact with the surface of the PMT window. Note that for 0j > sin'* (1 /Nq) = 42.4 deg, the 
photon will always be reflected. 

Assuming a uniform distribution over a hemisphere, a cap with a cone angle of 0 = 42.3 deg 
represents only 26 percent of the hemisphere’s surface area. Hence, 74 percent of the photons are totally 
reflected. In the modeling process, whenever the direction cosine of a photon measured with respect to 
the surface normal is 


do <d 


mm 


-y 


1 - (l/N G y 


(24) 


total reflection is obtained. 

If there is a top or bottom reflection, it is modeled by setting 

d 3 = -d 3 . (25) 
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PROBABILITY OF REFLECTION 
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Figure 5. Intensity of reflected light. 
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Figure 6. Geometry of reflection. 
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For a side reflection dj and d 2 are first transformed into components perpendicular and tangential to 
the side by 


dl' 


cos a 

sin a 


d l 

d 2 ' 


-sin ce 

cos a 


1 

<N 

1 


(26) 


with 


a = tan' 


1 


(u/v) 


The perpendicular component is then used for the total reflection comparison 




< d 


min 


If there is a total reflection it is achieved by reversing the sign of dj' and then transforming back to the 
unprimed system as follows: 


d l 


cos a 

-sin a 


-di' 

1 

a* 

ts> 

L 


sin a 

cos a 


I 

a* 

1 


Considering the endpoint as the startpoint of a new partial path, the calculations are repeated 
until one of two conditions occurs. 


1) The photon strikes the photocathode, liberates a photoelectron and the photon is absorbed. 
The photoelectron is counted if it strikes the effective photocathode area and a tally is kept. 

2) The photon strikes the sides at such an angle that it could never strike the photocathode, 
but would continue to be reflected indefinitely. The photon is ignored as soon as this condition is 
recognized. Given an effective photocathode radius of r and a window radius of R, the photocathode 
cannot be struck if (Fig. 7) 


|d 2 , l/d 1 ' > r 



r2 


(28) 
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Figure 7. Geometry of rejection. 


I. Modeling of the PMT Window Bottom 

When the PMT is manufactured, the photocathode material is sputtered uniformly over the 
bottom of the PMT window. The photocathode radius of 1.5 mm used in the model represents the 
effective area of the PMT electronics. Thus, all photons which strike the bottom of the window and are 
transmitted through liberate an electron from the photocathode material. Electrons produced in the 
effective area contribute to the photon count, while those outside this region make no contribution. 
Photons which liberate electrons are considered to be incapable of later producing another electron and 
are thus eliminated from the simulation after electron production. 


J. Photon Sink 

The 1 mm thick wall of the PMT which connects the window with the rest of the tube acts as a 
photon sink in that those photons which go down the walls are modeled as escaping from the system. 


K. Photon Absorption 

The intensity of light through 1.83 m of glass is reduced by 50 percent. Thus, the probability 
of a photon traversing through L mm of glass is 

p = e -kL 0 < L < oo (29) 

and 

P = 1 - kL for small L (0 < L < 1 9 mm) , 
where k = 0.000379 mm'l 


L. Photon Scattering 

Scattering of internal reflections was investigated, but was found to be negligible and was con- 
sequently ignored in the simulation. 
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ill. RESULTS 


Because of ST’s expected exposure to charged particles in its space environment, it is believed 
that the ST’s pointing capabilities would be degraded. This had been recognized as inevitable during the 
time in which the ST will traverse the charged particle dense South Atlantic Anomaly (SAA), and con- 
sequently, ST would be put temporarily on gyro-hold without FGS use. However, concern arose over 
the effect of the more energetic galactic cosmic ray nuclei that would be present in that part of the 
ST’s orbit outside of the Earth’s radiation belts. 

Since an attitude control simulation for the ST was under development, a logical step was to 
construct a model which simulated the photon noise caused by charged particles. This model, when 
installed in the ST attitude control simulation, would generate the photon counts in a random manner, 
whereby the impact of the charged particles on the allocated pointing error budget could be assessed. 

Based on previous work in the area of photon detection by Teich, et al. [15-18], the first 
approach was to investigate the applicability of the Neyman Type-A distribution and the Thomas distri- 
bution which belong to the category of cluster (self-exciting) point processes. Cluster processes are 
characterized by a primary (mother) process which generates at each point secondary (daughter) events. 
When only daughter events appear in the final process and when primary and secondary distributions are 
both Poisson, the resulting distribution is known as the Neyman Type-A counting distribution with 
probability density function 


P(n;p,s) 


OO 


£ 

m=0 


(s m) n 
n! 


e -sm 


p m e-P 
m! 


(30) 


The parameter p represents the rate at which primary events occur and s is the average number of 
secondary counts per primary. 

The Thomas distribution is similar, except that primary events are counted in the final process. 
The probability density function is given by 


P(n;p,s) 


Y (s m)n ~ m e -sm P me ~ P 
(n-m)! m! 

m=0 


(31) 


It was soon recognized that the particle’s path length through the PMT window [equations (17) 
and (19)] would be an important variable in any model attempting to estimate the photon counts caused 
by charged particles. As the setting became more complex, it was realized that a simulation would be 
required to unite all the relevant factors that would contribute to the photon noise. Secondly, the model 
could act as a forum from which the sensitivity of the parameters could be studied. 

A computer program was written (Appendix F) to simulate charged particles striking the PMT 
window and the resultant Cerenkov radiation and fluorescence. Stereo visual computer graphics were 
particularly helpful in checking both the mathematics and programming during the modeling process 
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(Appendix D). The stereo pair in Figure 8 depicts the PMT window as struck by a charged particle 
(proton in this case) with entry point A and exit point B. Conically directed lines in Figure 9 show the 
cone of luminescence comprising the Cerenkov radiation. Figure 10 shows the paths of the fluorescence 
photons which follow after the Cerenkov radiation with a time constant of about .1 ms. Note the internal 
reflection. The computer graphics also pointed out the case in which the photocathode was struck by the 
fluorescence photons but not the Cerenkov radiation, demonstrating that the pulse amplitude of the 
Cerenkov photon noise is not a reliable warning of the subsequent fluorescence noise. Furthermore, the 
converse situation was also observed. This points out one of the difficulties in applying the Neyman or 
Thomas distribution, since sometimes the primary event (Cerenkov noise) is present in the final count, 
and sometimes not. However, in most applications, this would not pose a serious problem, since these 
two density functions have very similar properties. 



Figure 8. Proton striking the PMT window. 



Figure 9. Cerenkov radiation caused by proton strike. 



Figure 10. Fluorescence caused by proton strike. 
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A. Photon Noise Caused by Cosmic Ray Protons and Alpha Particles 


Frequency statistics were computed for cosmic ray protons and alpha particles. Figure 1 1 shows 
the bivariate frequency distribution for Cerenkov radiation and fluorescence emitted in the PMT window 
caused by proton strikes. Figure 12a depicts the bivariate distribution for photon noise, i.e., those 
photons in Figure 11 which actually struck the effective photocathode region (Table 2 gives the 
numerical results). Figure 11, 12a, and 12b indicate poor correlation between Cerenkov and fluorescent 
yield at the photocathode and indicate poor efficiency for rejecting fluorescent events by the upper 
level threshold of the PMT PAD. Figure 12b shows the same information with the height magnified and 
the peaks at the origin clipped. These figures illustrate that very little of the radiation emitted along the 
charged particle’s path actually strike the effective photocathode. This can be attributed to the following: 

1) Small effective photocathode, i.e., it represents only one percent of the total area. 

2) Photon loss because of transmission, absorption, and photon sink. 

3) Photon loss because of the photocathode region outside of the effective area which make no 
contribution to the photon noise. 



Figure 11. Bivariate frequency distribution of Cerenkov radiation and 
fluorescence emitted along the proton’s path. 
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Figure 12a. Bivariate photon noise distribution. 



Figure 12b. Bivariate photon noise distribution magnified with peaks clipped. 
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TABLE 2. BIVARIATE FREQUENCY DISTRIBUTION OF PHOTON NOISE CAUSED BY 
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Figures 13a and 13b represent the probability density functions for Cerenkov noise and fluores- 
cence noise obtained for cosmic-ray proton strikes. The mean number of simulated Cerenkov counts 
caused by proton strikes is 0,37. For comparison with experimental results, the inset figure is included 
where those strikes which yielded zero counts were ignored (since the multiplier photo tube recognizes 
only those photons which strike the photocathode). Thus, the mean obtained when zeros are omitted 
is 2.75. The mean of the fluorescence density function is 0.423 (zeros included) and 2.21 when zeros 
are omitted. The numerical values of the PDF and CDF (probability and cumulative density function, 
respectively) are given in Table 3. 

Figures 14a and 14b are the simulated probability density functions for alpha particle photon 
noise, with means 1.47 and 6.34 (zeros omitted) for Cerenkov noise and 1.16 and 3.10 for fluorescence 
noise. Table 4 provides the numerical values of the noise PDF and CDF for alpha particle fluorescence. 
Figure 15 depicts the simulated probability density function for the length of the charged particle’s 
path through the PMT window. Figure 16 is the distribution of emitted Cerenkov photons caused by 
proton strikes and Fig. 17 is the distribution of emitted fluorescence. 

One proposed method for eliminating the effect of fluorescence is to recognize the Cerenkov 
pulse amplitude height (proportional to the number of Cerenkov photons which strike the photocathode) 
and then blank the PMT counter for the next millisecond, thereby eliminating the majority of the 
unwanted fluorescence noise. 

However, the simulation points out that it often is the case that fluorescence noise is produced 
without Cerenkov photons having struck the photocathode. This may be largely attributed to the directed 
nature of the Cerenkov radiation in which photons exit the window with little or no internal reflection 
but the fluorescence photons, in their more random flights, eventually strike the photocathode either 
directly or after reflection. 
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TABLE 3. FLUORESCENCE NOISE DISTRIBUTION 
CAUSED BY PROTONS 


PHOTON NOISE 

PDF 
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0 
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1 
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0.919515 

2 
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0.999366 

14 

0.000370 

0.999739 

15 

0.000261 

1.000000 


TABLE 4. FLUORESCENCE NOISE DISTRIBUTION 
CAUSED BY ALPHA PARTICLES 
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Figure 14a. Probability density function of Cerenkov noise caused by alpha particles. 
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Figure 14b. Probability density function of fluorescence noise caused by alpha particles. 
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Figure 17. Distribution of emitted fluorescence photons caused by proton strikes. 


B. SAA Protons and Alpha Particles 

The difference between protons and alpha particles trapped in the SAA and the cosmic ray nuclei 
is their energies. The cosmic ray nuclei’s energy is on the order of magnitude of a BeV whereas SAA 
particle energies are on the order of magnitude of a MeV. Particles with energy less than 300 MeV con- 
tribute no Cerenkov radiation [0 = 0 in equation (20)] while the stopping power term dE/dX appearing 
in equation (19) has its minimum mean value at around 1 BeV, increasing sharply for less energetic 
particles. 

Consequently, the amount of Cerenkov radiation will be quite small and fluorescence abundant 
during the time in which ST will traverse the SAA. 

As expected, simulation studies for both SAA protons and alpha particles resulted in negligible 
amounts of Cerenkov radiation and large amounts of fluorescence. 


C. Parameter Study of PMT Window Height 

Also of interest were the effects of changing the PMT window dimensions. Comparisons were 
made for cosmic ray protons with the PMT window radius remaining fixed and the height (thickness) 
varied. 


Since the number of Cerenkov photons and fluorescence photons is proportional to the charged 
particle’s path length through the window, it seemed intuitive that the number of photon noise counts 
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would increase with the increased window size. However, simulated means of Cerenkov and fluorescence 
monotonically decreased as the height increased. The apparent anomaly can be attributed to the fact 
that the window’s surface area increases more rapidly than the window’s linear dimensions, meaning that 
although more photons are produced by the charged particle, there is even less likelihood that they will 
strike the photocathode since the ratio of the area of the photocathode to the window surface area is 
decreasing at a rate faster than the mean charged particle path length is increasing. 

Since we are working only with the distribution of induced counts, conditional on having a 
charged particle strike, it is difficult to assess the value of this result since clearly the rate of charged 
particle strikes increases as the size of the PMT window increases. Thus, one would need to work out the 
multivariate problem when searching for an optimal PMT window size. 


D. Comparison with Experimental Results 

An experiment was conducted by the Space Science Laboratory at Marshall Space Flight Center 
to measure the photon noise with the use of a multiplier photo tube identical to the one we modeled. 
A beam of muons was selected by a separate detector in coincidence with the top of the PMT window 
and a mean of three counts was recorded for Cerenkov radiation and two for fluorescence. 

The appropriate parameters were adjusted in the model to simulate charged particles striking the 
top of the PMT window. Directions cosines of the particle’s path were modeled as follows: 

1) The angle 0 of incidence (angle measured from the window top normal) was assumed to 
follow a distribution 

f(0) = — cos 2 60 , 0 < 0 < 7t/12 

7T 

where experimental evidence suggested that f « cos*-^ [5], For mathematical convenience, we used the 
cos^ function. 

2) The azimuth angle was assumed to be uniform (0,27r). 

The simulated means were 4.5 and 1.7 for Cerenkov and fluorescence counts, respectively, zeros 
having been suppressed. The ratios of the simulated means for the space environment to the simulated 
means of the experimental set-up are 0.53 for Cerenkov counts and 1.33 for fluorescence. The experi- 
ments performed with muons in the laboratory have large estimated errors (±50 percent) in Cerenkov 
and fluorescence yields. Therefore, the means obtained from the simulation are compatible with the 
experimental values. An accelerator experiment using particles with a high charge number would be 
required to perform a definitive comparison. Such an experiment would yield an accurate value for the 
fluorescent efficiency. 

Simulated means for Cerenkov noise (C) and fluorescence noise (F) of the lab experiment and the 
space environment for cosmic-ray protons is given below: 
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Simulation of 

Lab Experiment Space Environment 



C 

F 

C 

F 

Zero counts 
included 

0.27 

0.23 

0.22 

0.43 

Zeros deleted 

4.5 

1.7 

2.4 

2.2 


The higher rate of Cerenkov radiation in the simulation of the lab experiment may be explained 
in that the charged particles are directed at the PMT window from above. Thus the Cerenkov radiation is 
directed conically towards the effective photocathode region. In the space environment, equal numbers of 
particles would be expected to enter the window through the bottom as well as the top. Cerenkov radia- 
tion emitted by charged particles entering from the bottom would be directed out of the top, with very 
little being reflected and then striking the photocathode. 


IV. CONCLUSIONS 


A stochastic model for Cerenkov radiation and fluorescence resulting from the direct interaction 
of a charged particle with the multiplier phototube window has been developed. This model simulates 
the bivariate probability density function for the number of Cerenkov photons and fluorescence photons 
which strike the effective photocathode region which contribute to photon noise. It is applicable to 
charged particles trapped in the Earth’s radiation belts and to the more energetic cosmic ray nuclei. 
The model’s parameters were adjusted to correspond with the physical environment of an experiment 
designed to estimate photon noise caused by cosmic ray protons. Experimental mean counts corresponded 
very well with the model’s prediction, thereby justifying a certain degree of confidence in the model s 

ability to predict. 

The density functions for fluorescence noise obtained from the simulation of cosmic ray nuclei 
and alpha particles (Tables 1 and 2) have been inserted into the ST Fine Guidance Sensor model to 
provide randomly generated photon noise in the FGS simulation. This allows for the evaluation of their 
effect on the overall ST pointing error budget. 

A simulation of carbon nuclei strikes was also conducted and the resultant bivariate histogram of 
Cerenkov noise and fluorescence noise clearly pointed out that the Cerenkov pulse amplitude height could 
not reliably be used to signal the associated fluorescence noise. A limit was therefore placed on the error 
of the Fine Guidance Sensor interferometer signal. This limit minimizes the effect of the strikes. 

Variation of the PMT geometry showed that the average fraction of the total number of photons which 
intercepted the photocathode decreased as the thickness of the PMT window was increased. This, together 
with an accurate model for the rate at which charged particles strike a PMT window as a function of 
window dimensions could be applied in the design of multiplier phototubes. The model will also be used 
to simulate various charged particle/glass interaction experiments in astrophysics. 
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APPENDIX A 


UNIFORM RANDOM NUMBER GENERATOR 


The advent of high-speed digital computers has made it possible to use simulation techniques 
incorporating probabilistic features. These simulations are generally called Monte Carlo simulations and 
are resorted to whenever the system being studied is not amenable to analytical methods. An integral 
part of these simulations is the use of random numbers having a certain specified distribution charac- 
teristic of the process being studied. 

The basic element of all Monte Carlo simulations is the uniform random-number generator. Once 
uniform random numbers are available, all other desired distributions can be obtained either by use of the 
probability integral transformation or by applying some known relationship between the desired distribu- 
tion to be generated and the uniform distribution. 

One of the most widely used uniform random number generators today is the congruential type 
first proposed by D. H. Lehmer in 1949 [19]. Accordingly, the random number generator takes the 
form 


X i+1 = (A X A + B) Mod M , i = 0,1,2,... 


(A-l) 


where A, B, and M are integers and X 0 is a “seed” number. Since the congruential relationship is 

cyclical, the sequence of random numbers will repeat after a certain period. The statistical behavior of 
the generated random numbers is predominantly governed by the choice of the multiplier A and the 
modulus M. Therefore, the most widely used generators are of the type B = 0. Many empirical and 
theoretical tests have been developed to assess the “goodness” of random number generators. Based on a 
test called the lattice test [20], Carroll [21] found the congruential generator with A = 29,903,947 and 

M = 2^ - 1 to be a satisfactory generator for the computer on which we would run the simulation. 
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APPENDIX B 


GENERATION OF RANDOM NUMBERS FROM NON-UNIFORM DISTRIBUTIONS 


Many of the generators for continuous distributions are a direct application of the probability 
integral transformation [22], For a given uniform random number U between zero and one, a random 
number X having a desired distribution F(X) can be obtained by solving the equation U = F(X) for X. 
Since this process requires the determination of the inverse of the cumulative distribution function, its 
practicality depends upon the availability of explicit expressions or convenient approximations for this 
inverse cumulative distribution function. Generation of random numbers from discrete distributions is 
handled in a manner analogous to that of continuous distributions. Given a uniform random number U 
and cumulative discrete distribution F(X), we seek the least value of X for which F(X) > U. This X is 
the desired random number having discrete probability density function f(x). 


32 



APPENDIX C 


UNIFORM DISTRIBUTION OVER THE SURFACE OF A SPHERE 


Lacking any information to the contrary, we have assumed that the direction of a charged 
particle’s path is uniformly distributed. This is sometimes called the postulate of equal apriori 
probabilities. The direction can be expressed by direction cosines and our goal was to select these in 
such a way that each surface element of a unit sphere has equal probability of containing the unit vector 
tip obtained from the direction cosines of the particle’s path. The azimuth angle </> is taken to be uni- 
formly distributed from 0 to 2n. The elevation angle is determined by observing that the area of a band 
(width determined by an angle 6 ) around the surface of a unit hemisphere is given by 2n sin 6 (Fig. C-l). 
The requirement that the surface elements be equi-probable gives 


Pr(0 < 6) = (2 7r sin 6)1 2 ir = sin 6 


(C-l) 


which is the cumulative distribution function for the elevation angle. From equation (C-l) we get 


6 = sin' 1 U 


(C-2) 


as the desired random elevation. To generate points on the surface of a sphere, we let U range from 
-1 to +1. Thus, the direction cosines are given by 


V 


\J 1 - sin^ 6 cos 0 

d 2 

= 

1 - sin^ 6 sin <p 

_ d 3_ 


sin 6 


where </> = 2 ttU] and sin 6 = 2(U2 - 0.5), where Uj and U 2 are independent random numbers uniformly 
distributed between -1 and 1. 

Figure C-2 demonstrates points uniformly distributed on the surface of a hemisphere. 
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APPENDIX D 


STEREO VISUALIZATION WITHOUT OPTICAL AIDS 
(Cross-eyed Stereo) 


On many occasions in engineering and physical analysis it is useful to be able to sketch in three 
dimensions. To fulfill this wish in many cases, we have used a convenient technique which requires no 
optical devices other than one’s eyes. All that is required is a stereoscopic pair of images. One additional 
capability is necessary. The observer must be able to cause the lines of sight of his eyes to converge; 
i.e., one must cross one’s eyes. The stereo projections are formed as shown in Figure D-l. The images 
are reversed and viewed as in Figure D-2. With a little practice, one can easily learn to reconstruct 
mentally, the 3-dimensional scene from the reversed stereo pairs. The interested reader should try several 
viewing distances (the farther away the page the less crossing of the eyes is required and the easier it 
becomes to focus the images). Squinting may also help as it increases one’s depth of focus. When one 
first looks at a stereo pair, one focuses on the page and sees two similar but separate images. As one 
begins to cross his eyes, the two images become four. Continue crossing the eyes until the interior pair of 
images come together. Since the line connecting corresponding points on the images must be at the same 
angle about the line of sight as the line connecting the eyes, it may be necessary to rotate the page or 
rock the head until these two images become superimposed and seem to merge into a stereo image. This 
technique does require some practice but once mastered it can be very useful for easy visualization in 
3 dimensions. 

If a computer with plot capability is available, we can construct the necessary stereo projections 
from a set of points and lines that represent the object of interest. We have referred to such a representa- 
tion as a wire frame model because of the appearance of the image. Let P be a representative point of 
the model. Each point P is projected into the picture plane S as shown in Figure D-3. The point P is 
projected to the eyepoint E and the line PE intersects the picture plane S at P'. P' is the projection of 
P onto S. The set of all points P' projected from object points P together with the connecting lines form 
the desired projection. We set up a reference frame in the plane S. To do this, we must specify which 
way is up (so to speak). Let _u u be a unit vector in this direction and u r = u u x is a unit vector in S 

pointing to the right. We place the origin of the S coordinate reference at O. Observe that j: Q = jg + d J; 

where r Q and _rg are position vectors of O and E respectively. From the geometry shown in Figure D-3, 
we can see that 


Ip' =IE + (ip - Ie) d /£ ■ (ip - Ie) 


From this we can compute 


xp' = (r p / - r 0 ) • u r 
yp' = (ip' - Iq) ‘ Hu 
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Since 


• u u = ' i^r = 0 , 


x p » = (r p ' - r E ) -_u r 


yp' = <IP' - £e) • «U 


The set of points (xp', ypO plotted conventionally forms the desired projection. Size can be altered by 

scale adjustments. These projections are then placed as desired. Also, the values used for d and eye 
separation s are arbitrary and can be adjusted for convenience or eye comfort. In real life s a 65 mm 
and d = 250 mm for comfortable reading; however, it may be more comfortable for d to be larger. 
Some initial experimentation with this technique should establish desireable settings. 



Figure D-l. Stereo projection. 
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RIGHT IMAGE 


LEFT IMAGE 


S 

Figure D-2. Stereo reconstruction. 



S 


Figure D-3. Projection geometry. 



APPENDIX E 


RESTRICTED ENERGY LOSS DISTRIBUTION 


The stopping power or mean energy lost per unit pathlength in a medium by a charged particle 
is proportional to the particle’s nuclear charge squared. However, the stopping power is only the average 
energy lost to the absorber. Identical particles with the same trajectory will produce a distribution of 
energy losses often referred to as “straggling.” This variation in energy loss is the result of the finite 
number of events that occur in the ionization energy loss process. 

The stochastic variation of the energy lost to a thin absorber by a charged particle was first 
attacked by Landau [23]. His derivation assumes that only a small fraction of the particle’s initial energy 
is lost and no upper limit is placed on the amount of energy which may be transferred to a single 
electron. 

Vavilov [24] obtained a more rigorous solution without the approximation of unbounded single 
energy transfers. 

Watts [25] makes the distinction that energy loss depends on the type of material and energy 
deposited depends on the surrounding geometry as well. 

Under the assumption that all electrons with energy below ej deposit their energy in the volume 

element where they are produced, energy loss becomes equivalent to energy deposition. Energies above 
ej are transported out of the region by scattering of electrons and are referred to as delta rays. 

The probability density function for energy loss A over a pathlength of x cm is given by [25] 
/c(l + ej / 3 2 7 ) 

— — — OO 

f(A ; x)=— — e max J e^ 8 ^^ cos [Xu + Kg 2 <u)] du (E-l) 

ne l 0 

where k = x|/ej, ej is the cutoff energy, and e max is the maximum possible energy transferred to an 
electron given by 



38 



where y is Euler’s constant 0.577215665; 


e l ft 

gi (u) = [In u - Ci(u)] - u Si(u) - cos (u) , (E-2) 

e max 

e l ft 

g2(u) = u [In u - Ci(u)] + Si(u) + sin u . (E-3) 

e max 

Other terms are as follows: 

Ci is the cosine integral approximated by: 



f y + In u 

u <0.01 



1 ' + taU+ 2 ; (2j,!2j 

0.01 < u < 1.0 


Ci(u) = < 


(E-4) 


G(u) sin u - H(u) cos u 

1.0 < u < 50.0 



sin u 
u 

u > 50.0 


and the sine integral is approximated by 




u 

u < 0.01 


Si(u) = < 

y, (_l>ju 2 j +1 

r-i (2j + l)!(2j + 1) 

J=0 

0.01 < u < 1.0 

(E-5) 


7t/2 - G(u) cos u - H(u) sin u 

1.0 <u< 50.0 



, cos u 
7r/2 

u > 50.0 



u 


where 

i /u^ + ai + a-?u 4 + a^u 2 + a A \ 

G(u) = - ( L- L - i ) (E -6) 

u V u° + a^u 6 + a^u 4 + a^u 2 + ag / 
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(E-7) 


U(ll) _ 1 / u8 + b l u6 + b 2 u4 + b 3 u2 + b 4 \ 


u 2 \ u^ + b$u b + bgu^ + byu 2 + bg / 


aj = 38.027264 

a 5 = 40.021433 

a 2 = 265.187033 

a 6 = 322.624911 

a 3 = 335.677320 

a 7 = 570.236280 

a 4 = 38.102495 

a 8 = 157.105423 

bj = 42.242855 

b 5 = 48.196927 

b 2 = 302.757865 

b 6 = 482.485984 

b 3 = 352.018498 

b 7 = 1114.978885 

b 4 = 21.821899 

b 8 = 449.690326 


The average stopping power A is given by the Bethe Equation as [27] : 


_ 4 it N n z ; 2 e^ Z p 

A = 

M e c 2 j3 2 A 


In 


2M e c 2 0 2 


*adj ^ P ) / l + 


2M„ 




© 


Z 2 


,2 \'P 


(E-8) 


where 


P = density of the material 

z = nuclear charge of the particle 

c = velocity of light in a vacuum 

e = electronic charge in electrostatic units (ESU) 

M e = rest mass of the electron 

A = the atomic weight of the stopping material 

N 0 = Avogadro’s Number 
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Z = atomic number of the stopping material 
I ac jj = the adjusted ionization potential 

P = ratio of the incident particle velocity to the velocity of light 
SjCj = the sum of the effects of shell corrections on stopping power 
5 = the polarization effect correction term 
Mp = proton rest mass 

The substitution of r 0 = e /M e C , where r Q is the classical electron radius, provides a useful 
change of variables in the evaluation of equation (E-8). For a particle of energy E/nucleon, (3 is given by: 


P = 



H — 

\1 + E/M p c 2 


) 


2 


Given the incident particle’s energy and length of trajectory through the glass, the PDF of the 
energy loss may be computed using equation (E-l). Thus, the stochastic variation of energy loss, or 
straggling, may be introduced into the Monte Carlo simulation by generating random energy losses from 
this distribution. 

However, these computations were too time prohibitive to evaluate for each simulated charged 
particle strike. It was thus decided to generate a family of these distributions as a function of charged 
particle energy, external to the photon noise simulation. Hence, the energy loss for each charged particle 
striking the PMT window could be quickly obtained from this family of distributions using numerical 
methods. 

The pertinent parameters in the energy loss distribution are pertinent to protons passing through 
silicon glass: 

ej = 1 MeV 

z = 1 

p = 2.33 gr/crn^ 

Z = 14 

I adj = 170 eV 

Mp = 1.67 x 10"^4 grams 


cutoff energy 
nuclear charge of a proton 
density of silicon 
atomic number of silicon 
ionizing potential of silicon [27] 
proton rest mass 
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electron rest mass 


Mg = 9.1 1 x 10' 28 grams 
N 0 = 6.02 x 10 23 
r = 2.82 x 10’^ 3 cm 


Avogadro’s number 
classical electron radius. 


The last two terms in the Bethe equation limit the stopping power for very high energies and may 
be set to zero in the present application because the cutoff energy ej effectively provides the same 

limiting function. For a detailed discussion of these terms, see Janni [27] and Sternheimer [28]. 

With these initial conditions set and A and /3 computed, the integrand in equation (E-l) may now 
be examined. 

The importance of making a careful examination of this surface is two-fold. First a suitable range 
for the energy loss must be determined so that the density function is zero for A < Aj and also for 

A > A 2 , where Aj and A 2 are two energies for which energies below Aj or above A 2 occur with 

negligible probabilities. Secondly, the upper limit of integration (the lower limit is always zero) must be 
chosen so that the PDF is smooth and that increasing the upper limit of integration provides no 
significant contribution to the integral. Unusual waves and negative values in the PDF are indications that 
this upper limit needs to be increased. 

The energy loss axis is constructed as 65 evenly spaced points between 0.5 and 4.5 (units are in 
MeV). These limits were chosen as discussed above and will vary with both charge particle type and 
absorber material. Thirty was found to be suitable as the upper limit of integration. The cutoff energy 
ej was determined to be 1 MeV [29], 

The integrand in equation (E-l) is illustrated in Figure E-l where the variable of integration u is 
in the u-direction and the energy loss A in the A-direction. Surfaces of this type proved to be invaluable 
in choosing the range of A and u for various particles and energies. 

Performing the numerical integration of the surface yields the desired PDF given in Figure E-2 
with its numerical values given in Table E-l. 

Vavilov’s distribution may be obtained by setting ej = e max in Equations (E-l) through (E-8). 
Figure E-3 compares Vavilov’s distribution with Watts’ restricted energy loss distribution. 

To implement the stochastic behavior of energy transfers into the photon noise model, a family 

of these distributions was evaluated for proton energies ranging from 700 MeV to 10 3 MeV. This is 
depicted as an energy surface in Figure E-4, where the E-axis represents proton energy (logjQ) and the 

A-axis represents deposited energy (logjQ), i.e., logjg (column labeled A in Table E-l). 

To implement this into the Monte Carlo simulation, a random proton energy E' is generated for 
each proton striking the PMT window. Referring to Figure E-4, the resultant PDF for energy transfer 
A, denoted by f(A;E'), is seen to be the intersection of the A-Z plane with the surface at E = E'. 

To obtain the desired energy transfer A, the equation 
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A 

f f (A;E') dA = U 
o 


is solved in terms of A, where U is a standard uniform random number. This gives A as an energy trans- 
fer randomly distributed according to Watts’ restricted energy loss distribution. 



Figure E-l. Integrand of restricted energy loss distribution [equation (E- 1 ) ] for 
1 BeV proton through 1 cm of silicon glass. 



ENERGY LOSS A (M»V) 

Figure E-2. Restricted energy loss distribution for 1 BeV proton 
through 1 cm silicon glass. 


43 



1 . 5 


UATTS 



UATTS 

MEAN SIGMA 
1.69 0.421 


VAVILOV 

MEAN SIGMA 
1.62 0.635 



VAVILOV 


Figure E-3. Vavilov’s distribution (set ej to e max = 3.3 MeV) versus Watts’ restricted energy 
loss distribution (ej = 1 MeV) for proton with incident energy of 
1 Bev through 1 cm of silicon glass. 
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Figure E-4. Family of energy loss PDF’s for protons of energy E through 

1 cm of silicon glass. 




TABLE E-l. 


A (MEV) 

£ (A; 1 cm) 

A (MEV) 

f (A, 1 ( 

0.5625 

0.0000 

2.5625 

0.1291 

0.6250 

0.0001 

2.6250 

0.1076 

0.6875 

0.0000 

2.6875 

0.0894 

0.7500 

0.0000 

2.7500 

0.0740 

0.8125 

0.0002 

2.8125 

0.0611 

0.8750 

0.0023 

2.8750 

0.0503 

0.9375 

0.0173 

2.9375 

0.0413 

1.000 

0.0737 

3.0000 

0.0339 

1.0625 

0.2030 

3.0625 

0.0277 

1.1250 

0.4081 

3.1250 

0.0226 

1.1875 

0.6521 

3.1875 

0.0183 

1.2500 

0.8802 

3.2500 

0.0148 

1.3125 

1.0508 

3.3125 

0.0119 

1.3750 

1.1474 

3.3750 

0.0096 

1.4375 

1.1744 

3.4375 

0.0077 

1.5000 

1.1469 

3.5000 

0.0062 

1.5625 

1.0834 

3.5625 

0.0049 

1.6250 

0.9995 

3.6250 

0.0039 

1.6875 

0.9069 

3.6875 

0.0031 

1.7500 

0.8137 

3.7500 

0.0025 

1.8125 

0.7249 

3.8125 

0.0019 

1.8750 

0.6427 

3.8750 

0.0015 

1.9375 

0.5684 

3.9375 

0.0012 

2.0000 

0.5020 

4.0000 

0.0009 

2.0625 

0.4426 

4.0625 

0.0007 

2.1250 

0.3888 

4.1250 

0.0006 

2.1875 

0.3398 

4.1875 

0.0005 

2.2500 

0.2948 

4.2500 

0.0003 

2.3125 

0.2536 

4.3125 

0.0003 

2.3750 

0.2164 

4.3750 

0.0002 

2.4375 

0.1834 

4.4375 

0.0002 

2.5000 

0.1543 

4.5000 

0.0001 
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APPENDIX F 


APL COMPUTER PROGRAM 


The following set of APL (a programming language) functions constitute the simulation model. 
The function HITS on p. 51 is the main program, calling all the others with NHT the number of hits 
desired. The function HTIC on p. 53 contains the values of most of the model parameters. Any time a 
function name ends with a A it contains a line-by-line commentary on the function of the same name 
without the A. Some functions have commentary within. 


V A<-S_ATAN C 

Cl) A+( 20Ci US*2'+C*2)*0 . \ S 

[2] n 

[3.1 n 4 QUADRANT ARC TANGENT 

V 

V Y+BC UiYl;Y2;YZ;UU 

[1] Il-EL+*(U<0.01 y >/U 

[21 Y2*-EL+(9UU)+SUM1 UU+( ( 0 . OlSff) a ( ffsiH /U 

[3] I3+((lOUU)*BG UU)-(20UU)*BH UU+( ( 1 <U) A ( U<5 0 ) )/U 

[4] ->0 -,Y+Y1,Y2,Y3, ( lOUU ) *UU+( U> 5 0 ) /U 
[51 n 

[6] fl EQ.EH APPROXIMATION TO COSINE INTEGRAL 

V 

V R+KE BETHE Z ; C ; MEV FERG ; A; ZM ; RHO ; IZ ; MPC S ; 5 ; K1 \K2 ; K3 


[1] ft KE-ENERGY/ NUCLEON OF CHARGED PARTICLE (MEV) 

[2] n Z-CHARGED PARTICLE CHARGE NUMBER 


[31 £N-2. 99792458#10n . . . .VELOCITY OF LIGHT IN V AC U U M - - CM / SEC 

[41 MEVFERG+il . 60219175 1 6r MEV FROM ERG'S 
[ 5 J MP+1 . 6 72 614ff_2 4fl . . . .PROTON REST MASS GRAMS 

[6] ME*-9. 28n ELECTRON REST MASS GRAMS 

[7] NO+6. 022169£23fl . . . .AVOGADROS NUMBER- - ATOMS / MOLE 


[81 R0+2.82E 1 3n CLASSICAL ELECTRON RADIUS-- CM 

[9] A *- 2 8.1fl ATOMIC WEIGHT OF SILICON 

[10] ZM*- 14n ATOMIC NUMBER OF SILICON 

[111 RHO*- 2.33fl DENSITY OF SILCON GRAMS /CM 

[12] IZ*-0. 00017fl IONIZING POTENTIAL OF SILICON-MEV 


[13] MPCS+MEVFERGxMP*CxCfi PROTON REST ENERGY - - MEV 

[14] MECS+MEVFERGxME*CxCfi ELECTRON REST ENERGY--MEV 

[15] B<-( l-( 1+KEiMPCS)* 2) *0 . 5 fl RATIO OF PARTICLE'S VEL . TO C 

[16] S*-1-B*B 

[17] Kl<-kxON0xMECSxR0xR0xZM*ZxZ 

[18] K2+( l+( (MEiMP)*2)+(2xME*MPxS*0. 511*0.5 

[19] K3+2xMECSxBxB*IZxS 

[20] R*-(KHBxBxA) x( ( -BxR) +9K3iK2 ) 

[ 21 ] n 

[22] ft THE BETHE EQUATION MEV/CM/GMxGM 

7 
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V Y+BG U i AA; NUM;DEN 

[1] AA+ 38.027264 265.187033 335.67732 38.102495 40.021433 
322.624911 570.23628 157.105423 

[2] NUM++/ (((pU) ,5)pl,4tj4i«)xJ7o.*8 6 4 2 0 

[3] DEli++/(. ( (pi/) , 5)pl, 4t44]x[/o.*8 6 4 2 0 

[4] +0',Y<-NUMiVEN*U 
[5J n 

[6] a £<3.£6 ffl(J/l 

V 

V 7^Bff U;BB;NUM-,DEN 

[1] BB-*-42 . 2 42855 302.757865 352.018498 21.821899 48.196927 
482.485984 1114.978885 449.690326 

[2] NUM<-+/( ( (pB) ,5)pl,4 + BB)xB°. *8 6 4 2 0 

[3] DEN++ /(((pB),5)pl, 41BB)xtfo.*8 6420 

[4] ->0 ;Y+-NUM*DENxUxU 

[5] A 

[6] A EQ .El G2(tn 

V 

V 7«-BS B;71;72;73;BB 

Cl] 71-<-( B< 0 . 01)/B 

[2] Y2+-SUM2 ( (0. 01<B''a(B<1') )/B 

[3] 73^(00. 5)-( (20UU''*BG UU'+(BH UU > x loi/B«-( ( 1<B> A ( B<5 0 U /U 

[4] ->-0 ; 7-«-7 1,72,73, (oo. 5 ^ - ( 20BB ) *UU+( U> 5 0 'I /U 

[5] A 

[6] a EQ.E 5 APPROXIMATION TO SINE INTEGRAL 

V 

V R+CDFI CDF 
[1.1 +(FLAG = 1 )/Ll 

[2] ->0 ; R++ / CDF<UNF 

[3] B1 : £-<-[ 0 . 5 +ST+ (ST*0 . 5 ) xNORM 

[4] A 

[5] A INVERSE CUMULATIVE DISTRIBUTION FUNCTION 

V 

V CR-,TR',D; DCR ; UVWCR ; XYZCR ; PR ; Z4P 

[1] CNTCR+ 1 x ( FS- 1 ' , FS=Q 

[2] CNTCRll + FSl + ({RSS>SPS XYZHT ) *XYZHTl 3 1 = 0 W ( RSS>SPS UVWHT 'I A 
UVWHT [31=0 

[31 +(B.B<l)/0 

[4] TR+-§ ( EUL2 20DHTL31 ) + . *EULZ DHTl2}ATAN DHTl 1] 

[51 LM:DCR+TR+ . x(EUL302*UNF' + . *(EUL1 2OiEN') + .*0 0 1 

[6] UVWCR+DCR EXITTB XYZCR<-XYZHT+LXYZHT*UNF 

[7] -*( DOPR-3 ) / BO 

[ 8] PR-t-PRC R<-PRC R ,UVWCP. , 0 ; PRC R+-PRC R ,XY ZC R , 1 

[ 9 ] BO ZAP= 1 ) /B1 -,CNTCRl 1 + FS1+CNTC Pi l+FSl + ZAP<r( RSS>SPS UVWCR'* A 
UVWCRi 3 1 < 0 . 01 

[10] CNTCRl 1+FSl+CNTCPll+FSl+RF UVWCR;D+DCR 
[111 DOPR- 0) / LI 

[12] PRCR+PR 

[13] LI :+( 0<NCR+NCR-1) / LM 

V 
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V CRL 

[1] n DATA IS SAVED SPLIT INTO TOP/ BOTTOM AND SIDES, 1 OCCUPIE 
S THE NOT USED SPOT 

[2] ft COUNT IF ENTRY OF EXIT POINT OF PARTICLE PATH IN CATHODE 

[3] n OMIT CERENKOV IF ENERGY TOO LOW 

[4] r TRANSFORMATION FROM COORD SYS WITH Z-AXIS ALONG MAIN PATH 

TO XYZ COORD SYS, USING CERENKOV PHASE ANGLE (' 20*BN'i OF A 
BOUT 47 DEG 

[5] r GET DIR COS WITH CERENKOV ANGLE AND UN IF AZIMUTH ABOUT MA 
IN PATH 

[6] a PICK UN IF DISTR ALONG PATH AND GET EXIT POINT 

[73 A SKIP NEXT LINE FOR NO GRAPHICS 

[83 r SAVE GRAPHICS DATA 

[93 A ZAP= 1 IF CATHODE HIT, ENDS PATH OF THIS PHOTON BY JUMP TO 

LI 

[103 R GET COUNT FROM REFLECTIONS AND SAVE 
[113 ft SKIP NEXT LINE IF NO GRAPHICS 
[123 ft SAVE GRAPHICS DATA 

[13 3 ft REPEAT FROM LM UNTIL NCR= 0 (NO OF REMAINING CERENKOV PHOT 
ONS ^ 

V 

V R+EUL 1 A 

[13 R+ 3 3 p ( 1 0, lO-A,-A, A+O0. 5 3 [ 1 2225423 53 
[23 ft 

[33 ft EULER ANGLE ROTATION MATRIX ABOUT X-AXIS 

V 

V R+EUL 2 A 

[13 R* 3 3 p ( 1 0, 10-4,-4,4+00. 5)[5 2321242 53 

[23 ft 

[33 ft EULER ANGLE ROTATION MATRIX ABOUT Y-AXIS 

V 


V R+EUL 3 A 

[13 R + 3 3 p ( 1 0,10-4,-4,4+00. 5>[5 4235222 1] 

[ 2 j ft 

[3] ft EULER ANGLE ROTATION MATRIX ABOUT Z-AXIS 

V 

V UVW+H EXITSD XYZ ; X ; Y ; Z ; U ; V ; AZ ; SC ; D2 

[ 1 3 X+XYZ [13; Y+XYZ [23; Z+XYZ [ 3 3 ; D2+ 1 DHT [ 2 3 

[23 U* 2 xRMxDHT [ 1 3 x D2 i 1 -DHT [ 3 3 x DHT [ 3 3 

[33 V+(D2 *UiDHTtl~\)-RM 

[43 LZ+DHTl3lxUiDHTlll 

[53 +(A Z<HM-Z) /NTS 

[63 +TB ;SC*-(HM-Z) iLZ;H+HM 

[73 NT S :+( LZ>- Z) / NBS 

[8] SC+-ZihZ-,H+0 

[93 TB:U+SCxU‘,V+-(SC*V + RM''-RM 

[103 +0 -,UVW+U,V,H 

[113 NBS-.UVW+U ,V ,Z+ AZ 

[12] ft 

[13] ft EXIT POINT FOR SIDE ENTERIES 

V 
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V UVW+D EXITTB XYZ;A;B;C;SC;AXYZ 

[1] H^-(I7Z[3]-HMxZ)[3]>0)tD[3] 

[2] UVW*-XYZ + AXYZ*-D*H 

[3] +(RMS>SPS UVW )/ 0 

[4] A+SPS A XYZ 

[5] B+-+/ (XYZxLXYZ)ll 2] 

[6] SC+(B+( (B*B) + (RMS-SPS XYZ 1 * A > * 0 . 5 1 * >3 

[7] UVW+XYZ + AXYZ+SCx&XYZ 

[ 8 ] Ft 

[9] n EXIT POINT FOP TOP/ BOTTOM ENTRY 

V 

V FL ;D;PRi DFL ; ZAP 

[1] CNTFL+ 1*(FS=1\FS=0 

[2] -+■ ( 0 -NFL ) / 0 

[3] PT :XYZFL+XYZHT+AXYZHT*UNF 

[4] DFL+GENDC 

[5] UVWFL+DFL EXITTB XYZFL 

[6.1 CNTFLL 1+FSl+CNTFLl l + FSl + ZAP^-(RSS>SPS UVWFL'i MJVWFLl 3 ]<0 . 01 

[7] ->( DOPR-O ) /L0 

[ 8 ] PR+PRFL+PEFL ,UVWFL , 0 ;PRFL+PRFL , XYZFL , 1 
[91 LOiD+DFLi+(ZAP=l)/Ll 

[10] CNTFLL l+FSl+CNTFLl 1+FSl+RF UVWFL 

[11] +(DOPR=0)/Ll 

[12] PRFL+PR 

[13] LI : + ( 0<NFL+NFL-1 ^ / PT 

V 

V FLA 

[1] n SEE COMCRL 1 ] 

[2] n JUMP OUT FOR NO FL PHOTONS 

[3] a PICK UN IF DISTR ALONG PATH 

[4] a GET UN IF DISTR DIR COS 

[5] a GET EXIT POINT 

[6] a DETERMINES IF CATHODE IS HIT AND SAVES COUNT 

[7] R SKIP NEXT LINE FOR NO GRAPHICS 

[8] Ft SAVE GRAPHICS DATA 

[9] Ft TAKE CARE OF NEXT PHOTON IF CATHODE IS HIT 

[10] Ft SAVE FLUORESCENT COUNT 

[11] Ft SKIP NEXT LINE FOR NO GRAPHICS 

[12] Ft SAVE GRAPHICS DATA 

[13] Ft COUNT NLF DOWN BY ONE AND REPEAT FROM PT IF NOT ZERO 

V 

V D+GENDC ; AZ ; S 

[1] L^-( (( 1-5x51*0. 51 x( 20/^Zl , 10^Z^-02x[/]7F'' ,S^2x 0.5 + UNF 

[ 2 J Ft 

[3] Ft RANDOMLY DISTRIBUTED DIRECTION COSINES 

V 
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7 GENKINEN 
[IT +(SAA=1)/L 

[ 2] ENPN+JFBEV*KEBEV*-( 1-UNF*( l-CUTOFF*l-GAMl > '>*il-GAMl 
[31 +LL 

[4] L:ENPN<-JFMEV*GENLOGNOPM 4.7 In MODEL FOR SAA PROTONS AND A 
LPHA PARTICLES 

[51 LL : BN*-N7 056*(1-(1 +ENPN i M0C2 ) * 2^ *0.5 

C 6 1 a 

[7] fl KINETIC ENERGY AND BETA ANGLE 

V 

V E+GENLOGNORM PARAM ;MU iV AR 

[11 VAR+PARAMl 2! iMU+PARAMl 11 
[2] E+*MU+(VAR*0. S)*NORM 

[31 A 

[41 fl GENERATES N LOG-NORMAL RANDOM NUMBERS WITH PARAMETERS 
[51 fl MU AND VAR ( MU AND VAR ARE THE MEAN AND VARIANCE OF 
[6] fl THE DEFINING NORMAL DISTRIBUTION) . 

V 

7 GENRTRPAV \N\TH ; CTH ; STHT ; CTHT \RT ;RP ;TT ;TP 
[11 N*-N 7056 

[2] CTH+DMIN+(l-DMIN)*0,0.02x\50 

[31 TH+ 20 CTH 

[4] CTHT*- ( l-STHT*STHT*-NxlOTH)*0 . 5 

[51 RT+RTxRT*-(.CTHT-N*CTH) ZCTHT+NxCTHfi REFLECTION , PARALLEL 

[61 RP+RpxRP+(CTH-NxCTHT) iCTH+NxCTHTfi REFLECTION , PEPPENDICULA 

R 

[71 TT*-(CTHTt N xCTH ) xTTxTT*-2xNxCTHiCTHT+NxCTHpi TRANSMISSION , PPP 

[81 TP+(CTHTiNxCTH) xTPxTP*-2xNxCTHiCTH+NxCTHTn TRANSMISSION , PAP 

ALLEL 

[91 pppp^y^i + o. 5xpr+pp fi of /?r+pp 

v 

V Y^Gl U 

[1] ->-0 ;!-*-((( pB) , pP") p (-( 200) +0x55 £/ ) ) + ( RxBxEl i E' « . x ( ®U ) - BC U 

[21 n 

[31 fl EQ.E2 

7 

7 Y*-G2 U 

[11 +0 iY*-( (BxBxEUE) ° . xBS U) + ( ( pP ) , pU ) p ( lOP )+Ux ( ®P ) -BC U 

[21 fl 

[31 fl EQ.E 3 

7 
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V HITS NHT ; DHT ; XYZHT ; UVWHT ; NCR ; NFL 

[1] HTIC 

[2] PRCNTHT+- \ 0 

[3] PRLD+PRCNTCR+PRCFTFL + 2 OpO 

[4] LhPRHT+PRCR+PRFL+H OpO 

[5] HT 

[6] CR 

[7] FL 

[ 8 ] PRCNTCR+PRCNTCR , CNTCR-, PRCNTFL+PPCNTFL ,CNTFL 
[9] + (0<NHT+NHT-D / LI 

V 

V HITS A 

[1] ft INITIALIZATION 

[2] ft tfP 

[3] n DATA SAVING 

[4] n RESET GRAPHICS DATA 

[5? ft GENERATE RANDOM PARTICLE STRIKE PATH 

[6] ft EVALUATE CERENKOV PHOTONS 

[7] ft EVALUATE FLUORESCENT PHOTONS 

[8] ft SAVE CENRENKOV AND FLUORESCENT COUNTS 

[9] ft GO BACK FOR NEXT HIT UNTIL HITCOUNT NHT IS ZERO 

V 

V HT ; TR •, L ; ST 

[1] DHT+GENDC 

[2] ■+(RTB>UNF' 1 / TB 

[3] XYZHT+O A-RM) % RM*UVF\FS+1 

[4] UVWHT+HM EXITSD XYZHT 

[51 PRCNTHT+PRCNTHT ,(RSS>SPS UVWHT ) aUVWHTZ 3 1 < 0 . 01 

[61 TR+EUL302*UNF 

[71 XYZHT+TR+ .*XYZHT 

[81 UVWHT+TR+. xUVWHT 

[91 +CL 

[ 101 TB:XYZHT+(RM*2 x 0 . 5+UNF ,UNF) >HM*DHT [ 3 1 < 0 ; FS+0 
[111 +(,RMS<SPS XYZHT) /TB 
[121 UVWHT+DHT EXITTB XYZHT 

[131 PRCNTHT+PRCNTHT , ( (RSS>SPS UVWHT) *UVWHTl3 1<0 . 01 1v (RSS>SPS 
XYZHT) AXYZtfr[3l<0. 01 

[141 CL: !+■(+/ UXYZHT<-UVWHT-XYZHT)*2)*0 . 5 
[151 PRLD+PRLD ,L,DHTl 3] 

[161 GENKINEN 
[171 -*■( BN^l ) / L0 

[181 NCR+CDFI POISSON ST+-KANC Rx-L*l -tBN*BN 

[191 ft L0 :NFL<-CDFI POISSON ST+KANFL* L*BN BETHE Z 

[201 L0 : NFL+CDFI POISSON ST+KANFL* L* WATTSPR ENPN*MEVFJ 

[211 -*( DOPR= 0 4 / 0 ; ft ' NCR ,NFL- ' ;NCR ,NFL 

[221 PRHT+-PRHT , UVWHT , 0 ; PRHT+PRHT , XYZHT . 1 

V 
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V HT A 

[1] a p^prippp p^r# pjp cop, PffrC3]>o is upward hit 

[2] ft GO TO TB IF TOP /BOTTOM HIT 

[3] A GET RIGHT OF SIDE ENTRY, SET SIDE FLAG 

[4] ft GET EXIT POINT FOR SIDE ENTRY 

[5] A SAVE MAIN STRIKE, 1 FOP PHOTOCATHODE HIT 

[6] a TRANSFORMATION ABOUT Z AXIS WITH RANDOM ANGLE 

[7] A TRANSFORM ENTRY POINT 

[8] A TRANSFORM EXIT POINT 

[9] A GO TO CL 

C101 n PICK TOP ( DHT C 3 ] > 0 1 OR BOTTOM POINT IN A 2 *RM SQUARE ; RES 

ET SIDE FLAG 

[11] A REDO POINT IF OUTSIDE RM RADIUS CIRCLE 

[12] A GET EXIT POINT 

[13] A SAVE A 1 IF ENTRY OR EXIT POINT HIT PHOTOCATHODE , OTHERWI 
SE SAVE 0 

[14] a GET PATHLENGTH 

[15] A SAVE PATHLENGHT AND DHT [ 3] ( UPHIT DHT l 3]>0> 

[16] A GENRE ATE KINETIC ENERGY 

[17] a SKIP CERENKOV PHOTON NCR GENERATION IF ENERGY TOO LOW 

[18] a GET CERENKOV PHOTON NO WITH POISSON AVERAGE 

[19] A GET FLUORESCENT PHOTON NO WITH POISSON AVERAGE ( FOR LOW Z 
) 

[20] A GET FLUORESCENT PHOTON NO WITH POISSON AVERAGE ( FOR HIGH 

Z) 

[21] A SKIP NEXT LINE IF NO GRAPHICS 

[22] A SAVE GRAPHICS DATA 

V 


52 



V HTIC 

[1] Z«-lf» PARTICLE CHARGE HUMBER 

[2] DOPR+lfi DORR - 1 FOR GRAPHICS ( OVERFLOWS WORKSPACE CAPACITY F 
OR LARGE Z) 

[3] SAA+On SAA= 1 RESULTS IN SAA CHARGED PARTICLES ; 0 IMPLIES COS 
MIC RAYS 

[ 4 1 n ******** DEDX2 NEEDED FOR SAA = 1, CTM/l/GP 1/1/ VtfTV NEEDED ** 

***** 

[5] il/7 0 5 6«-l . 4 8 7 fl PPPff I/l/PET 

[61 RHO+ 2.3 3iq DENSITY OF GLASS GRAMS /ML 

[71 i4*-2 8 . 0 8 6 fl ATOMIC WEIGHT OF SILICON 

[8] ZM+14fl ATOMIC NUMBER OF SILICON 

[ 9 1 JZ-*1 7 0n IONIZING POTENTIAL OF SILICON-- MEV 

[101 P0-*2 . 81 7 945’ 13rt CLASSICAL ELECTRON RADIUS- -CM 

[ill /1/0-6 . 02486E23R AVOGAGROS NUMBER ATOMS/MOLE 

[121 MP+- 1. 67261£' 24n PROTON REST MASS- - GRAMS 

[13] ME+9 .10956E 28fl ELECTRON REST MASS- -GRAMS 

[14] MCSQ+Q . 51097a ONE RYDBERG 

[151 DMIN+- ( 1 - ( *fl/7 0 5 6 ) * 2 ^ * 0 . 5fl COS OF TOTAL REFLECTION ANGLE 

[ 16] DMINM+-DMIN ; KRT+50 * 1 - DM IN ; GENRTRPAV 

[17] n SID+RFD * Ofl SIGMA OF REFLECTION VARIATION 

[18] KTR+O . 000379n TRANSMISSION LOSS, ABSORPTION 

[19] LO:+LO+Z 

[20] -+L1 ; GAMl-*-2 . 7 n Z = 1 PROTON GAM1 IS INTEGRAL FLUX OF PROTON 
INTENSITIES 

[21] +L1 ;GAMl+2. 7n Z = 2 4LPP/4 NASA PP-243, P 0.135 

[22] +Ll;GAMl+2 . Z=3 LITHIUM 

[23] -+L1 ; GAM 1+2 . 5 6r Z = 4 BERYLLIUM 

[24] +LliGAMl+2.3Sft Z = 5 B0B0J7 AS27P0 . JOURNAL , 2 1 8 : 2 6 9 - 2 8 5 , 1977 
NOV 15 

[25] +L1 -,GAMl<-2 . 76iq Z = 6 CARBON 

[26 ] ->-Ll ;0yiAfl«-3. 0 7fi Z = 7 NITROGEN 

[27] -+L1 ;G4Ml<-2 . 67iq Z = 8 OXYGEN 

[28] -*£1 ;G,4Ml*-3 . Oln Z = 9 FLUORINE 

[29] -*-Ll;0j«Jfl«-2.88fl_Z = lO 0£0ff 

[30] £1 :M0G2*-1. 673E 2 7x299792458*2 

[31] CVTOFF+lQOn CUTOFF FOR HIGH ENERGIES-- UNITS IN BEV 

[32] PM*- 9 . 5 n PM2 7 WINDOW RADIUS IN MM 

[33] RS+1.5H PHOTO CATHODE RADIUS IN MM 

[34] HM+ 1.5fl PMT WINDOW HEIGHT IN MM 

[35] RR+RM - In PIWG RADIUS IN MM, 1 4S MEASURED BY US 

[36] KLAM+ 244n 7056 WAVELENGTH INTEGRAL 

[37] RTB+RM* RM+HM 

[ 38] RMS+RMxRM’,RSS+RSxRS’,RRS+-RRxRR 

[39] KANCR+- (ZxZxo2t137. 03602'' xKLAMn FIRST FACTOR IN EQU . FOR CER 
ENKOV RADIATION (TN D-81471 

[40] KANFL+2 x 8 n FACTOR OF TWO IMPLIES 4 PI STERAD IANS (PERS . CO 
NVER. WITH WALTER VIEHMANN 1 

[41] r FACTOR OF 8 FROM TABLE 1 ( TRIALKALI 70 56'' OF TN P-8147 
[421 C IRC T+RSx CIRC 

[43] JFBEV+1 . 60219P lOfl JOULE FROM BILLION ELECTRON VOLTS 

[44] JFMEV+JFBEV * 1000 

[45] MEVFJ+iJFMEV 
7 
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V o+v INTEGRATE T 

[1] 0*- (F[2]-K[l])x(-o. 5 x ( (14-p2 , ''p,2 , ') + (14p2'')p , <j>[ 1 iT-f-tyrl + + /T+T 

[2] A 

[31 n INTEGRATES A MATRIX OR VECTOR WHERE T=F(V> 

V 

V Y+XY INTP0L1 XiNX;NY;Nl;N2iXl;X2’,YUY2 

[11 -*•( (n[l;l]>L/lW(r/I)>"lUJ[l;lVS 

[21 N2+1+N1++/ (<St(NY ,NX^ pX)> ( (NX+pX) ,NY+( pXY' [ 2 1 * pZY[ 1 ; 1 

[3] Z1+-JY[ 1 ;tfl] ;Z2-<-ZY[ 1 ;ff2l ;Y1«-JY[2 ;/fll ;72^7[2; W21 

[4] Y«-Yl + (Y2-Yl)x(Z-Zl WZ2-Z1 

[5] ->-0 

[61 E:'X OUTSIDE RANGE' 

[71 n LINEAR INTERPOLATION 

[81 fl ZY[ 1 ; 1 MUST INCREASE MONOTONIC ALLY , 

[9.1 fi BUT DOES NOT HAVE TO BE EVENLY SPACED 

V 


V R+NORM 

[11 R<-(( ' 2>CUNF)*0'.5)x2OO2xUNF 

[21 n 

[31 n NORMAL DISTRIBUTION BASED ON THE POLAR METHOD 

V 

V CDF+-POISSON ST;P0;P;K 

[l! -*-(Sr>17 7)/Ll 

[21 FLAG+0 

[31 CDF<-P0+*-ST 

[41 K+-P + 1 

[51 L0 : P*-Px ST *K 

[61 CDF+CDF, ( HCDF)+P*P 0 

[71 Jf-f-Z+l 

[8] ->( 0 . 9999> liCDF)/L0 

[ 9 1 -►O 

[101 LI :FLAG*-1 ;CDF*-\ 0 

[111 fl 

[121 n GENERATES POISSON DISTRIBUTION WITH MEAN ST 

V 


54 



Cl] 

[ 2 ] 

[3] 

[4] 
C5] 
[ 6 ] 
C7] 
[ 8 ] 

[9] 

CIO] 
Cll] 
C 1 2 ] 
C 13] 
C 14 ] 
C 15 ] 
C 16 1 
C 17 D 
C 18] 
C 19 ] 
C 20 ] 
C 21 ] 
C 2 2 ] 
C 2 3 ] 
C 24 ] 
C 2 5 ] 
C 26 ] 
C 27 ] 
C 2 8 ] 
C 29 ] 
C 30 ] 
C 3 1 ] 
C 3 2 ] 
C 3 3 ] 


Cl] 

C 2 ] 

C 3 ] 

C 4 ] 

C 5 ] 

C 6 ] 

C 7 ] 

C 8 ] 

C 9 ] 
CIO] 

Cll! 
C 12 ] 
C 1 3 ] 
C 14] 
C 15 ] 
C 1 6 ] 
C 17] 
C 18 1 


NN PROFAMILY UE ; LM ;E ;P ',C ; D ; A ;KE ; K;T ;U ; R 
KEV+PROTONKEV ; N+NN C 1 ] ; M+NN C 2 ] ; A>1 

V+B PTS 0.6 Hil+UZm-UllhU+M PTS IE 1 0 ,UE ; R+( N , pKEV ) p 0 

L : LM+KE BETHE Z«-l ;I«-pXF«- 2+(2 *K'>+KEViMM+ 2t ( 2*K) + i p REV 

T-*-o2*N0xMECSxR0 xR0xZxZtB*B 

E+2xMECS*BxB* ( 1~B*B ' xl + ( 2*S) +S*S+MEiMP 

P+TiEl+UEL+O .5772156649 

C-*-( iOEl)x*Pxl+ElxB*BxELiE 

D-*- ( ( ( 7 o . xjpl )-(N ,I)pLMU UlM>pffl >-(N ,1 'l pPx (B*B^+1 + (®F1 iE > - 
EL 

i?C ',MMl + { (N , I) pC) xINTi (N ,1 ,M) p(*(P° . *Mpl^xGl U ) ) x 20 ( D» . *10 + ( 
N ,1 ,M) p(P° . *Mp 1 ^ X(J2 U 
-*•( ( ( pFF7>*2 ]>F«-F+1 V£j^ 

PR0T0NFAMILY+R* ( ( 1+pff ) pi ) • . xy INTEGRATE^ 

•►0 

ft A? IS THE NO. OF PARTITIONS OF THE VARIABLE U OF 
ft INTEGRATION; N IS THE NO OF FREQ. PTS. FOR EACH 
ft OF THE RESULTANT DISTRIBUTIONS ; PP0T0NKEV IS THE 
ft VECTOR ENERGIES PER NUCLEON OF THE CHARGED PARTICLE ; 
ft R IS THE MATRIX WHICH CONTAINS THE FAMILY OF DENSITY 
ft FUNCTIONS ; UE IS THE UPPER LIMIT OF INTEGRATION ; 
fl A IS THE INTEGRATION STEP SIZE ; 5 15 OF 

ft THE PROTON'S VELOCITY TO LIGHT , I .E . , B = ViC; KEV IS 
ft PRIMARY PARTICLE ENERGY IN MEV /NUCLEON ; Z IS PRIMARY 
ft CHARGE NO.', NO IS AVOGADRO'S NO.', R0 IS CLASSILAL 
ft ELECTRON RADIUS i ME IS ELECTRON REST MASS', MP IS 
ft PROTON MASS (Z = l>; E IS MAX. POSSIBLE ENERGY 
ft TRANSFER', EL IS EULER'S CONSTANT ; T IS ZETA 
ft DEFINED IN EQ. FI; FI IS DELTA RAY CUTOFF ENERGY. 

ft 

ft 

ft PROFAMILY PRODUCES A FAMILY OF RESTRICTED ENERGY LOSS 
ft DISTRIBUTIONS . EACH GENERATED DISTRIBUTION CORRES- 
ft PONDS TO AN ENERGY DEFINED EXTERN ALLY BY THE VECTOR 
ft PROTONKEV WHICH SHOULD CONTAIN AN EVEN NUMBER OF 
ft ELEMENTS. 


CNT+RF UV W ', XY Z ', Z AP ’, N C 0 
CNT+NC0 + 0 

TBSD :->-(( UVWl 3 ] *0 ) *UVW[ 3 1 *HM) / SD 

+ (.DMIN> |PC3] )/TR* TOTAL REFLECTION 

+ (RTRPAV C T KRT*DMINM+ | Dl 3] 3*UNF ) /On TRANSMISSION 

TR :+L0 ;D*-1 1 1 *D 

SD:NC0<-NC0 + 1 

AL+UVW121ATAN UVW Cl] 

D+iEULZ AL)+.*D 
+ (DMIN<D CUV Oft TRANSMISSION 

•*( ( I Dl 2] )>£>C 1 ?*RSi (RMS-RSS) *0 . 5 'S / 0 ft SIDE REFLECTION MISSES 

PHOTOCATHODE 

D+(EUL3-AD + . xFx"l 1 1 

LO-.UVW+D EXITTB XYZ+UVW 

CNT+CNT+ZAP+(RSS>SPS UVW^ hUVWl 3 ]<0 . 01 

-►( (F7VC3]<0. 01)vJIZC3]<0. 0 1 ) / 0 ft NO BOTTOM REFLECTIONS 
+ (DOPR = 0 VF1 
PR+-PR ,UVW , 0 

FI :'•*■( UNF> 1 -KTRx ( + / ( UVW-XYZ ) * 2 ") * 0 . 5 V 0ft TRANSMISSION LOSS 
-»•( ( Z4P= 0 )AFC0< 10 )/TBSD 



V RF A 

[13 A SET COURT AND MAX SIDE STRIKES TO ZERO 

[23 A IF NOT ON TOP OR BOTTOM , GO TO GO 

[33 A IF I02MI REFLECTION , GO TO TR 

[43 A JUMP OUT IF TRANSMISSION (UN IF PICK FROM TRANSMISSION CUR 

VE) 

[53 A RFLECT AND GO TO 10 

[63 A INCR SIDE STRIKE NO 

[73 A GET SIDE INCIDENCE ANGLE 

[83 A TRANSFORM D WITH SIDE INC ID ANGLE 

[93 A JUMP OUT IF TRANSMISSION 

[103 A JUMP OUT IF SIDE REFLECTION MISSES CATHODE 

[113 A REFLECT AND TRANSFORM BACK TO XYZ COORD SYS 

[ 123 A SET PREVIOUS 'EXIT' POINT TO 'ENTRANCE' POINT AND GET NEW 
'EXIT' POINT 

[133 A CHECK FOR CATHODE HIT AND SAVE COUNT 
[143 A JUMP OUT IF NO BOTTOM REFLECTIONS 
[153 A SKIP NEXT LINE IF NO GRAPHICS 
[163 A SAVE GRAPHICS DATA 

[173 A JUMP OUT IF TRANSMISSION LOSS ALONG PHOTON PATH 

[183 A REPEAT FROM TBSD IF NO STRIKE AND LESS THAN 10 SIDE REFLE 

CTIONS 

V 

V R+AXIS SMOOTH S 

[13 F-+or s 

[23 F<-Ft( ( ltpff) pi) o . xAXIS INTEGRATE R 
[33 A 

[43 A SMOOTHES FAMILY OF VAVILOV DISTRIBUTIONS 

V 

V R+SPS I 

[13 i?W(I*lUl 23 

V 

V Y+SUM 1 U iNUMiDENiJ 

[13 NUM+(U° . *2xj) x ( ( pU) , pi) p" 1*2 | 2 5 

[23 DEN+UpU) ,pJ)p('.2*J)x2*J 
[33 +0 ;Y++ /NUMiDEN 

[43 A 

[53 A INFINITE SERIES WHICH APPEARS IN THE COSINE INTEGRAL 
[63 A APPROXIMATION IN EQ . F4 

V 

V Y+-SUM2 U ; NUM DEN ; J 

[13 NUM-*- (U ° . *l + 2xl)x((o0) , pi) p" 1*2 | J+~ l + i 25 
[23 DEN-*-( ( p U) , pi) p (.' 1 + 2*1) x 1 + 2 xj 

[33 -+0 ; Y+-+ / NUMi DEN 

[43 A 

[53 A INFINITE SERIES WHICH APPEARS IN THE SINE INTEGRAL 
[63 A APROXIMATION IN EQ.E5 

V 
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V U+UEF 

Cl] U+1E 9 x 1+?1000000001 

C2] n 

[3] r STANDARD UNIFORM RANDOM NUMBERS 

V 


V NN V AV PROTON UE ;T \E ;P ;C ; D ; S ; KE ;U iPROTONSUR ; LM 

[1] I*-pKE-*-10 0 0 ; N+-NNL 1 ] ; M+NNl 2 ] fl ENERGY IS IN MEV 

[2] LM+PP BETHE Z«-l 

[3] 7«-P PrS 0.5 4.5;P«-M PI’S 1 1P'10,Z/P 

[4] y -<-0 2 xPO *MEC S* R0 xR0xZ*ZiBxB 

[5] P«-2 xMECSxB*B* ( 1 ~B* B) x 1 + ( 2 *S) +S*S+ME*MP 

[6] P<-T*P1«-1 ;PZ>0 . 5772156649 

[7] C+( iOEl ) x*Px1+E1xBxBxELtE 

[ 8 ] ZM ( (7°. xJpl)-(jy,J)pLWH(iy,J]pffl>-(y,J]pPx(pxflWl + («ffl*g'- 

EL 

[9] PROTON SUR+( (N , I ,M) pC ° . 1 ) x ( ( N , J , p ( * ( p o . xMp D *G1 !/')x 20 ( 

D° .*U)+(N ,1 ,M)p(P° .xMpl)xG2 U 

[10] ft 

[11] R PRODUCES THE FUNCTION P(J?,A'> BEFORE INTEGRATION 

[ 12] fl 

[13] PDF+V SMOOTH U INTEGRATE PROTON SUR 

V 

7 DEDX+W ATT SPR EiU;I;J-,Z 

[1] I-*- 10 1 2 + + /P^I<P 

[2] U+UNF;J<-1;Z+- lO 

[3] R+2 ,pDEAXPR 

[4] L: Z+Z, (RpIVAVPR [;![</]] ,DEAXPR) INTP0L1 U 

[5] 4> e 7-<-e7'+l ) / L 

[6] DEDX+0 . IxRHOx ( 2 4pP4J[ J] , Z ' INTPOL 1 Pfi MEV/GP2/CM TO MEV / MM 
[71 fl 

[8] n ENERGY DEPOSITION DEDX GIVEN PARTICLE ENERGY E 

V 
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